Question 1

sf_uscities = readr::read_csv("D:/LifeInUCSB/Study/GEOG176A/week2/assignment5/geog-176A-labs/data/uscities.csv") %>% 
  filter(city ==("Palo")) %>% 
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
  st_transform(5070) %>% 
  st_buffer(5000) %>% 
  st_bbox() %>% 
  st_as_sfc() %>% 
  st_as_sf()

Question 2-3

st = list.files("D:/LifeInUCSB/Study/GEOG176A/week2/assignment5/geog-176A-labs/data", full.names = TRUE, pattern = "TIF")

s = stack(st) %>% 
  setNames(c(paste0("band", 1:6)))


cropper = sf_uscities %>% 
  st_transform(crs(s))

r = crop(s, cropper)

#Step 3

#The dimensions is 7811 rows and 7681 columns, 6 layers. 
#The CRS is WGS84.
#The cell resolution is x=30, y=30 (meters). 

#Step 4
#The dimensions is 340 rows and 346 columns, 6 layers. 
#The CRS is WGS84.
#The cell resolution is x=30, y=30 (meters). 
#R-G-B (natural color)
par(mfrow = c(1,2))
plotRGB(r, r = 4, g = 3, b = 2,stretch = "lin")
plotRGB(r, r = 4, g = 3, b = 2,stretch = "hist")

#NIR-R-G (fa) (color infared)
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 4, b = 3,stretch = "lin")
plotRGB(r, r = 5, g = 4, b = 3,stretch = "hist")

#NIR-SWIR1-R (false color water focus)
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 6, b = 4,stretch = "lin")
plotRGB(r, r = 5, g = 6, b = 4,stretch = "hist")

#My choice
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 7, b = 1,stretch = "lin")
plotRGB(r, r = 5, g = 7, b = 1,stretch = "hist")

#Colored stretch is a method to process different feature in map to keep a visual track. This method adjusts the brightness of a certain wave band to improve its identification and avoid confusion with other types of data.

Question 4

Q4.1

ndvi = (r$band5- r$band4) / (r$band5 + r$band4)

ndwi = (r$band3 - r$band5) / (r$band3 + r$band5)

mndwi = (r$band3 - r$band6) / (r$band3 + r$band6)

wri = (r$band3 + r$band4) / (r$band5 + r$band6)

swi = 1 / (sqrt(r$band2 - r$band6))

stack = stack(ndvi, ndwi, mndwi, wri, swi) %>% 
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
palette = colorRampPalette(c("blue","white","red"))

plot(stack, col = palette(256))

#description
#These pictures show the same distinction between floods near river and land, but with different coloration.
#Flood area in NDWI,MNDWI and WRI plots are depicted in red. In NDVI and SWI plots, floods are drawed in blue.

Q4.2

thresholding1 = function(x){ifelse(x <= 0, 1, NA)}
thresholding2 = function(x){ifelse(x >= 0, 1, NA)}
thresholding3 = function(x){ifelse(x >= 1, 1, NA)}
thresholding4 = function(x){ifelse(x <= 5, 1, NA)}

flood1 = calc(ndvi,thresholding1)
flood2 = calc(ndwi, thresholding2)
flood3 = calc(mndwi, thresholding2)
flood4 = calc(wri, thresholding3)
flood5 = calc(swi, thresholding4)


floodstack = stack(flood1, flood2, flood3, flood4, flood5) %>%
  setNames(c("NDVI-Cells less than zero", "NDWI-Cells greater than zero", "MNDWI-Cells greater than zero", "WRI-Cells greater than 1", "SWI-Cells less than 5"))
plot(floodstack, col = "blue")

Question 5

set.seed(09072020)

value = getValues(r)

dim(value)
## [1] 117640      6
value = na.omit(value)

#Q5.2
#The dataset has 117640 rows and 6 columns.

scvalue = scale(value)


kmeans12 = kmeans(scvalue, centers = 12, iter.max = 100)
kmeans_raster = r$band1
values(kmeans_raster) = kmeans12$cluster
plot(kmeans_raster)

table = table(values(flood1), values(kmeans_raster))

which.max(table)
## [1] 6
func_kmeans = function(x){ifelse(x == 3, 1, 0)}
cal_kmeans = calc(kmeans_raster, func_kmeans)
floodstack = addLayer(floodstack, cal_kmeans)
plot(floodstack, col = palette(256))

Question 6

6.1

kabletable = cellStats(floodstack, sum)
knitr::kable(kabletable, caption = "total area of the flooded cells", col.names = ("Number"))
total area of the flooded cells
Number
NDVI.Cells.less.than.zero 6666
NDWI.Cells.greater.than.zero 7212
MNDWI.Cells.greater.than.zero 11939
WRI.Cells.greater.than.1 8469
SWI.Cells.less.than.5 15201
layer 28924